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ABSTRACT 

Wc propose a new ray-tracing algorithm to measure the weak lensing shear and con- 
vergence fields directly from TV-body simulations. We calculate the deflection of the 
light rays lensed by the 3-D mass density field or gravitational potential along the 
line of sight on a grid-by-grid basis, rather than using the projected 2-D lens planes. 
Our algorithm uses simple analytic formulae instead of numerical integrations in the 
computation of the projected density field along the line of sight, and so is computa- 
tionally efficient, accurate and straightforward to implement. This will prove valuable 
in the interpretation of data from the next generation of surveys that will image many 
thousands of square degrees of sky. 
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1 INTRODUCTION 



Weak gravitational lensing (WL) is a promising tool to 
map the matter distribution in the Universe and constrain 
cosmological models, using the statistical quantities pri- 
marily constructed out of the observed correlations in the 
distorted images of distant source galaxies. In 2000, four 
teams announced the first observational detections of cosmic 
shear (Bacon et al 2000; Kaiser, Wilson & Luppino 2000; 
van Waerbeke et al. 2000; Wittman et al. 2000; Maoli et al. 
2001). Since then improved observational results have 
been published (Hoekstra et al. 2006; Fu et al. 2008; 
Schrabback et al. 2010), and it has been extensively used 
to investigate key cosmological parameters such as the 
matter density parameter f2 m , and the normalisation of 
the matter power spectrum as as well as for constrain- 
ing neutrino mass (Tereno et al. 2009). Much theoreti- 
cal progress has also been made in assessing the utility 
of cosmic shear in, for example, estimating the equation 
of state of dark energy w (Bridle & King 2007; Li et al. 
2009; Crittenden, Pogosian & Zhao 2009), as well as its 
role in testing theories of modified gravity (Schmidt 2008; 
Zhao et al. 2009, 2010a, b; Song et al. 2010) and constrain- 
ing quintessence dark energy (Chongchitnan & King 2010). 

On linear scales, one can use linear perturbation theory 
to calculate the WL observables for a given cosmology, such 



* E-mail: b.li@damtp.cam.ac.uk 
f E-mail: ljk@ast.cam.ac.uk 
X E-mail: gong-bo.zhao@port.ac.uk 
Sj E-mail: hz4@st-andrcws.ac.uk 



as the shear power spectrum or the aperture mass statis- 
tic, and compare these predictions to observational data to 
constrain the model parameters. However, the observables 
on nonlinear scales, which cannot be predicted theoretically 
without the help of TV-body simulations, can also provide 
valuable information to prove, or falsify cosmological mod- 
els. Making such predictions using TV-body simulations be- 
comes increasingly important as we move into a new era 
in weak lensing using large observational surveys. The next 
generation of cosmic shear surveys, e.g., the Dark Energy 
Survey (DES; www.darkenergysurvey.org) will be more than 
an order of magnitude larger in area than any survey to date, 
covering thousands of square degrees, and using several fil- 
ters that alfow photometric redshift estimates for the source 
galaxies to be derived. These surveys have the potential to 
map dark matter in 3-D at unprecedented precision, testing 
our structure formation paradigm and cosmological model. 

To obtain the statistics for WL from the outputs of TV- 
body simulations, one needs to construct numerous virtual 
light rays propagating from the source to the observer. By 
tracing these light rays along the lines of sight (l.o.s.), one 
could in principle calculate how much the original source 
image is distorted, and magnified. 

Conventional ray-tracing algorithms generally project 
the matter distribution along the paths of light rays onto 
a series of lens-planes, and use the discrete lensing approx- 
imation to compute the total deflection of the light rays 
on their way to the observer (Jain, Seljak & White 2000; 
Hilbert et al. 2009). The lens planes could be set up either 
by handling the simulation outputs after the TV-body simu- 
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lation is completed or by recording corresponding light cones 
on-the-fly (Heinamaki et al. 2005) and projecting later. Al- 
though this algorithm is the most frequently used in the 
literature, it requires a large amount of data, such as par- 
ticle positions, to be stored, and this would be difficult for 
simulations with very high mass resolution or very big box 
sizes, which are increasingly more common today. Further- 
more, projecting particles onto a number (~ 20 — 30) of lens 
planes will inevitably erase the detailed matter distribution 
along the lines of sight and oversimplify the time evolution 
of the large scale structure. 

One can also perform the lensing computation during 
the TV-body simulation process to obtain the projected (sur- 
face) density and/or convergence field directly (White & Hu 
2000). This method avoids the expensive storage of dump 
data at numerous redshifts and allows the detailed matter 
distribution to be probed. However, it does involve numer- 
ical integrations in the calculation of the projected density 
field and therefore certain overheads, because in order to 
make the integrals accurate one has to sample the density 
field rather densely. 

Motivated by the promise of cosmic shear surveys, and 
the need to make predictions of observables on nonlinear 
scales using cosmological simulations, in this work we in- 
troduce a new algorithm to preform ray-tracing on the fly, 
which is based on that of White & Hu (2000). We calculate 
the deflection of a light ray as it goes through the iV-body 
simulation grids using the 3-D density field inside the grids, 
instead of using the density field projected onto discrete 2- 
D lensing planes. Furthermore, the numerical integration 
is replaced by some exact analytic formulae, which could 
greatly simplify the computation. We will show our result 
in comparison with the fitting formula, and discuss how our 
algorithm can be applied to particle or potential outputs 
recorded in large simulations, and how we can go beyond 
the Born approximation and include the lens-lens coupling 
effect. 

This paper is organised as follows. We will introduce 
our algorithm in the next section, describe our simulation 
and present the results in Sect. 3, and close with a section of 
discussion and conclusion. Although we do not include lens- 
lens coupling and corrections to the Born approximation in 
our simulations, we will outline in Appendix A how these 
can be done. For simplicity, we shall consider a spatially 
flat universe throughout this work, but the generalisation to 
non-flat geometries is straightforward. We shall use "grid" 
and "grid cell" interchangeably to stand for the smallest unit 
of the mesh in the particle-mesh iV-body simulations. 



2 METHODOLOGY 

In this section, we will first briefly review the traditional 
'plane-by-plane' ray-tracing algorithm, and then detail our 
improved 'grid-by-grid' prescription. 



2.1 Conventional Ray-tracing Algorithm 

We work in the weak-lensing regime, meaning that the light 
rays can be well approximated as straight lines (Mellier 1999; 
Bartelmann & Schneider 2001). The metric element is given 



by 

ds 2 = a [(1 + 2$)dr 2 - (1 - 2$)dx ■ dx] (1) 

where a is the scale factor normalised so that a — 1 today, 
t is the conformal time, $ is the gravitational potential and 
x the comoving coordinate. We use units such that c = 1. 

Then the change of the photon's angular direction as it 
propagates back in time is (Lewis & Challinor 2009) 

f(x.) " 6 = "2 1^ Zz^Vpdx (2) 

Jo xxs " 

in which \ is the comoving angular diameter distance, £ is 
the angular position perpendicular to the l.o.s., £o = £(x = 
0) , Vj- denotes the covariant derivative on the sphere with 

respect to £ and <E> = <E> [x, the gravitational potential 
along the l.o.s.. The 2x2 distortion matrix is given by Aij = 
Vi^j = V£ oi £j(x), where £ch is the i-th component of £o, and 
is equal to 

Ay ee -2 r V foi V e . $ ( X , e) d X + 5 %J 

Jo X 

« -2 r ^%r^ v «- v ^ $ (* d * + 5i > W 

Jo X 

with i,j = 1, 2 running over the two components of £, and 
g{ X ,Xs)= {Xs ~ x)x . (4) 

Xs 

Note that to obtain Eq. (3) we have made the approximation 
V^ oi w , which means that lens-lens coupling is ignored. 
We shall discuss how to go beyond this approximation in 
Appendix A. 

This matrix is related to the convergence k and shear 
components 71,2 by 

. _ ( 1 - K - 71 -72 - U \ , v 

A -{ -72+a, l-K + 71 J (b) 

where w stands for the rotation, and 7 = Oyf + 7! ) 1 ^ 2 the 
shear magnitude. In the weak-lensing approximation, once 
the convergence is obtained, the shear is determined as well, 
therefore in practice we only need to compute k, 

k = 1~{A\+A\) /2= g( X ,Xs)V&dx. (6) 
Jo 

Under the Limber approximation (White & Hu 2000), the 
two-dimensional Laplacian in Eq. (6) can be replaced with 
the three-dimensional Laplacian, because the component of 
the latter parallel to the l.o.s. is negligible on small angular 
scales (Jain, Seljak & White 2000). Then, using the Poisson 
equation 

v 2 $ = \n m H 2 - (7) 

2 a 

where 8 is the matter overdensity and a the scale factor, we 
can rewrite Eq. (6) as 

k(£) = ^l m H% £ S g( X , Xs) S -^d X (8) 
in which we have written explicitly the a;x-dependence of k 
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(the ^-dependence is integrated out in the projecting pro- 
cess). Eq. (8) is the starting point of most ray-tracing sim- 
ulations. 

The most commonly-used ray-tracing method is the dis- 
crete lensing approximation. In this approach, the density 
field is projected onto a number of lensing planes (usually 
~ 20 — 30), and the light rays are treated as if they were 
deflected only by these plane lenses. Correspondingly, the 
term g(x,Xs) m Eq. (8) is evaluated only at the positions of 
these planes. 

The method of White & Hu (2000) incorporates the in- 
tegration in Eq. (8) directly into their iV-body simulation 
code, and performs the integral at every time-step. To re- 
alise this, N\os straight l.o.s. are generated to be traced. The 
rays have specified origin (the observer at redshift 0), open- 
ing (e.g., 3° x 3°) and orientation. As the iV-body simulation 
process evolves to the source redshift z 3 , the convergence is 
computed along each line of sight using Eq. (6) or Eq. (8). 
The l.o.s. integration is then carried out numerically for each 
time step, during which the photon travels from Xi to Xf> 
where the subscripts i and / literally stand for initial and 
final respectively, and hence they are used as the integration 
boundaries. The integrand gS/^Q/x 2 in Eq. (6) or g5(x) in 
Eq. (8) is considered to be constant during each time-step, 
and the integral is approximated by summing over all the 
time steps. The time sampling has to be sufficiently fine so 
as to guarantee the required numerical accuracy. 

One advantage of this algorithm is that k is computed 
step-by-step on the fly, so one can avoid the expensive disk 
storage required for storing particle dumps and the time- 
consuming postprocessing analysis. Moreover, in this ap- 
proach, there is no difficulty to make ultra- fine time sam- 
pling - the number of time slices can be as many as the 
number of time steps for the simulation (after z s ), which 
is a mission impossible for the postprocessing approach - 
making the result more accurate than the postprocessing 
approach. 

However, one does have to carry out the numerical inte- 
gration in Eq. (6) or Eq. (8), and to make the result accurate 
one has to sample the value of the integrand very densely 
(e.g., ~ 100 sampling points are dynamically chosen for each 
time step), which might cause certain overheads when a large 
number of light rays are traced and ultra-fine time-stepping 
is used. 



2.2 Improved Ray-tracing Algorithm 

In this work, we propose an improved ray-tracing algorithm 
by computing the convergence, shear and projected density 
fields on the very grid cells on which the iV-body simulation 
is performed. In our grid-by-grid approach, the l.o.s integra- 
tion can be carried out analytically, making the computation 
more efficient and accurate. Also, the light rays are deflected 
by the detailed matter distribution exactly as seen in the N- 
body simulations, making the ray-tracing and iV-body sim- 
ulations consistent with each other. A detailed derivation of 
the relevant formulae is given in Sections 2.2.1 & 2.2.2, and 
the basic idea is as follows. Take Eq. (8) as an example, the 
integrand is gS(x). Since our particle-mesh (PM) code auto- 
matically computes <5(x) on the regular mesh, the value of 
<5(x) at any point can be obtained by interpolation, and in 
particular we can compute the value along the line of sight 



as a function of the comoving distance x> an d the values of 
<5(x) at the vertices of the grid containing the said point . 
Note that the vertices themselves are regular grid points, 
and the values of <5(x) on the vertices are known. Using cer- 
tain interpolation schemes, trilinear, for example, 5(x) can 
be expressed as a polynomial of x, thus the integral can be 
carried out analytically. Therefore, no numerical integration 
is needed to compute n. Similarly, our algorithm can also 
be used to compute the integral of Eq. (G) analytically, as 
detailed in Section 2.2.2. 

Note that when the algorithm is applied to the density 
field <5(x), there is some subtlety, and this will be clarified 
in Sect. 2.2.1. 



2.2.1 Method A 

To compute the projected density field, we need to integrate 
along the l.o.s., and in practice this integration can be car- 
ried out progressively along the segments of lines of sight 
within individual cubic grid cells. The reason for such a pre- 
scription will become clear soon. Throughout this subsection 
we will use z to denote the coordinate rather than redshift. 

Fig. 1 shows two examples of such configurations, in 
which the part AB of a line of sight lies in a grid (see the 
figure caption for more information). The density value at a 
given point on AB could be computed using trilinear inter- 
polation, as long as we know the values at the vertices, which 
are denoted by p xyz (x, y, z = 0, 1). To be more explicit, let 
us define 
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(9) 

Suppose the point on AB we are considering has coordinate 
(x, y, z), then the density value is given by 



p(x,y,z) = c + ciAx + c 2 Ay + c 3 Az + c^AxAy 
+c 5 AyAz + c 6 AxAz + c 7 AxAyAz, 



(10) 



where 
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(11) 



where L denotes the size of the cubic cell, (xo, yo,zo) is the 



1 Because the line of sight is by approximation a straight line, 
once the comoving distance X to a point is known, the correspond- 
ing x,y, z-coordinates of that point can be expressed in terms of 
X and orientation angles, which are fixed when the l.o.s. are as- 
sumed to be straight. 
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Figure 1. Two examples of the linc-of-sight (lines) crossing a cubic cell of the simulation box at points A and B. The segment AB lies 
inside the cube. P X yz (x, y, z = 0, 1) are the eight vertexes of the cube. Projecting AB into the plane perpendicular to 2-direction and 
passing point A, then ip is the angle between the projection and x-direction, and is the angle between AB and that plane. For the 
given line-of-sight and cube, A, B, 9, ip are known or can be computed easily, and we also know the density values at the eight vertexes; 
we then want to integrate the density field along AB (or part of it). 



coordinate of vertex Pooo, Xa is the \ value at point A, and 
a, b, c are the coordinates of point A relative to Pooo • 

Because we express p(x, y, z) in terms of \ only, the line 
integral along AB could be rewritten as an integral over x 5 
and we have 
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(12) 



in which xi ^ Xa and Xu ^ Xb~ are the lower and upper 
limit of the integral respectively, x = X~XA, Xu = Xu~XA, 



2 Note that A and B axe the intersections between the l.o.s. and 
the grid cell, and not necessarily the two ends of the l.o.s. in one 
time step. But the integration is carried out for each time step, 
and so we do not always have \l = XA an d Xu = Xb ■ 



Xi = Xi ~ Xa and we have also defined 

di = c + ~ (aci + bc 2 + cc 3 ) 
Li 

+ (a6c4 + bcc 5 + acce) + -r^abccT, 
d'z = — cos 9 cos ipci + — cos 9 sin ipC2 + — sin 9c A 

Li Li Lj 

+ — cos 9 sin ip (ac4 + ccs ) 
L 2 

+ -pr COS 9 COS lb (feC4 + CC6 ) 

L 2 

+ sin 9 (bc 5 + ace) + -i- sin 9abcT 
L 2 

+ — cos 9 sin tpacc7 + — cos 9 cos ipbccr, 

1 2 1 

d$ = — cos 9 sin ip cos tpC4 + — sin 9 cos 9 sin ipc^ 
L 2 L 2 

+ -K; sin 9 cos 9 cos tpce H — ^ sin 9 cos 9 sin ijjaci 
L 2 L 6 

+ -rrs sin 9 cos 9 cos ipbcr 
L s 

1 2 
+ ~znr cos 9 sin lb cos ibcci , 
L s 

d,4 = -j-rr sin 9 cos 2 9 sin tb cos ifjc-j. (13) 

Li 

Note that, by writing the result in the above form, we have 
separated the treatments for four types of variables: 

(i) a, b, c, xa, Xu.i'- a, b, c, xa are determined by the direc- 
tion of the light ray and the specific grid cell under consid- 
eration, and Xu,i depend only on the considered time step 
and xa- Note that a,b,c must be determined carefully, and 
for each grid cell at least one of them vanishes, but exactly 
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which of them vanishes varies from ray to ray and from grid 
cell to grid cell; 

(ii) 9, ip: these specify the direction of the light ray, and 
terms involving them only need to be computed once, i.e., 
at the beginning of the simulation, for a given line of sight; 

(iii) C0-7 - these are determined by the values of p at the 
vertices of a grid, and must be evaluated for each grid that 
the light ray passes through; 

(iv) L,\ s : these are constants for a given simulation. 

Therefore once p xyz is known, the integral can be performed 
analytically without much computational effort. This is not 
unexpected, because once the density is known at the ver- 
tices of the grid, we should know the density at any point 
inside the grid using interpolation, and no more information 
is needed to carry out the integral. If we consider a different 
grid, a different set of p xyz needs to be used, and this is why 
our algorithm is based on the individual grids. 

There are two technical points which need to be noted. 
First, in Eq. (12) p should be replaced by p/a in practice. 
It is true that a could be expressed as a function of \ as 
well once the background cosmology is specified, but this 
will lead to more complicated expressions. Therefore in our 
simulations we simply take a to be constant during each 
time step. This is certainly only an approximation, but we 
should note that a is considered as constant during each 
time step in the iV-body simulations anyway. Indeed, as we 
see in Section 2.2.2, the factor - does not appear if we use 
V|5> instead of p in the integral'. 

Second, as has been mentioned by various papers (e.g. 
Jain, Seljak & White (2000); White & Hu (2000)), the use 
of the three dimensional Laplacian [Eq. (8)] instead of the 
two dimensional one [Eq. (6)] is at best an approximation. 
We have to test the validity of this approximation. In fact, 
as we show below, the error caused by this approximation is 
actually not negligible. To see this, recall that 

k = r g (V 2 - V 2 ) <£>d X 
Jo 

= ^n m H 2 Q £ S g 5 -d x - [gV x $>]f 

+ I gV x $d x + I g'V x $d X (14) 
Jo Jo 

in which a prime (overdot) denotes the \ (time) deriva- 
tive, and the last three terms come from the treatment of 
V x $, including integration by parts. The common argu- 
ment is that the second term actually vanishes as g = 
and V x $ < oo at \ — Xs an d X = 0, and the last two 
terms are negligible. This is true in the ideal case, but while 
our algorithm [and that of White & Hu (2000)] is applied 
the second term is no longer zero because of the following 
reasons: 

(i) It is unrealistic to make the simulation boxes big 
enough to contain the whole light cone, and in practice peo- 
ple tile different simulations to form a complete light cone. 

3 This just reflects the fact that during each time step of the N- 
body simulation, the 1/a factor in the Poisson equation is treated 
as constant. The nature of numerical simulation (discreteness in 
time) dictates that we cannot do better save decreasing the length 
of time-steps, which we cannot always keep doing in reality. 



Unless a periodic tiling of the same box is adopted, we ex- 
pect the matter distribution and thus the potential $ to 
be discontinuous at the tiling boundaries. As a result the 
second term in Eq. (i i) should read 

[sV x <£>]* s = lgV x ^ n + [gV x <f>]™ +■■■ + [gV x $]™ 

in which Xu,i correspond to the values of x when the light ray 
goes through a given box, which is labelled as 1, 2, • • • , N. If 
the matter distribution is smooth at the boundaries of the 
boxes, then V x $(x = Xll) — V x <&(x = Xu2) and so on, so 
all terms cancel. However, if the matter distribution is not 
smooth, as is the case for many tiling treatments, then such 
cancelling will not happen and [pV^]*" will turn out to be 
nonzero in the numerical calculation although it should be 
zero in theory. 

(ii) Using the same argument as above, we could find that 
this discontinuity problem appears not only on the bound- 
aries of the tiled simulation boxes, but also at each time-step 
in the simulations and each time when the light ray passes 
through a grid of the simulation box. For the former case, 
suppose that during one time step the l.o.s. ends at point C, 
then C is also the point where this l.o.s. starts during the 
next time step. However, the values of V x $ at point C are 
generally different in the two time steps because particles 
have been advanced, and so a discontinuity appears. For the 
latter case, our piecewise l.o.s. integral and the interpolation 
scheme dictate that the values of V x $ at a point D on the 
interface of two neighbouring grids could depend on which 
grid is supposed to contain point D (remember the inter- 
polation scheme uses the values of V x $ at the vertices of 
the containing cell), and naturally a discontinuity in V x $ 
appears at the interface of the two grids. Note that these 
discontinuities are inevitable due to the nature of numeri- 
cal simulation (the discreteness in time), and decreasing the 
grid size or the length of time steps does not help because 
then such discontinuities will only appear more frequently 1 . 

The way to tackle these problems is as follows: we know 
that [<7V x $]q s vanishes rigorously in principle but is nonzero 
because of the nature of the simulation; meanwhile, the same 
discontinuity problem also appears when calculating the first 
quantity on the right-hand side of Eq. (14). The errors in the 
numerical values for these two quantities are caused by the 
same discontinuity and could cancel each other. The exact 
value of this error can be obtained by computing [pV x $], be- 
cause this quantity is zero in theory and its nonzero value is 
completely the error. In our simulations, we compute [<?V X $] 
explicitly whenever the light ray passes a grid, and subtract 
it according to Eq. (14): this way we can eliminate the error 
in the integration of gS/a due to the discontinuities. 

As for the third and fourth terms in Eq. (14), the third 
term is nonzero but small in reality, but in our simulations 
it vanishes because <3? is assumed to be constant during any 
given time-step. This will cause certain unavoidable errors, 
that we anyway expect to be small. The fourth term has 
as small a contribution, but fortunately we can perform the 

4 Interestingly, the discrete lens-plane approximation does not 
have this problem (as long as simulation boxes are tiled periodi- 
cally so that matter distribution is smooth on the tiling bound- 
aries), because it does not treat the l.o.s. integral on a grid- by-grid 
basis. 
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integral exactly and analytically as we have done for the first 
term in Eq. (14). 

We have run several tests to check the accuracy of the 
approximations, and found the following: 

(i) If we simply replace the two dimensional Laplacian in 
Eq. (6) with a three dimensional one, as in Eq. (8), then the 
difference is of order 10% and even much larger for the rays 
for which \k\ is small. Note that Eq. (6) can be evaluated 
exactly as will be described in Section 2.2.2. 

(ii) If we explicitly calculate the term [gV x $]p s for each 
cell crossed by a ray, and subtract it according to Eq. (14), 
the difference between Eqs. (6, 8) is brought down to the 
level of 1-2%. 

(iii) If we further include the contribution from the fourth 
term of Eq. (14), the difference will fall well within the per- 
cent level. 



2.2.2 Method B 

The method described in Section 2.2.1 is only applicable to 
Eq. (8), while there are also motivations for us to consider 
Eq. (6). For example, the use of the three-dimensional Lapla- 
cian instead of the two-dimensional Laplacian in Eq. (8) is at 
best an approximation and only works well on small angular 
scales. This is even worse in the discrete lensing approxima- 
tion, because the photons of equal distance from the ob- 
server are certainly not in a plane but on a spherical shell, 
and this has motivated more accurate treatments such as 
the prescription proposed by Vale & White (2003). As an- 
other example, within the current framework the shear is 
not computed directly but from its relation with ft. There is 
certainly no problem with this, but it will be even better if 
we can compute 71,2 directly and compare with the results 
obtained from k. 

Our generalised treatment here is quite simple, taking 
advantage of the fact that the particle-mesh codes also give 
us the values of 5>(x) and (if necessary) ViVj$ at the regular 
grid points. For simplicity, let us assume that (1) the central 
line of sight is parallel to the x-axis, and (2) the opening of 
the lines-of-sight bundle is a square with its sides parallel to 
y, z-axes respectively. In the two-dimensional plane perpen- 
dicular to the line of sight, the i = 1,2 directions are set to 
be longitude and latitude respectively. We also define 
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= V m V x $, 
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EE VyVyQ, 
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EE V,V Z $, 
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= VxVyi" = 






= V y V,$ = 




■uj 


EE V,V 2 $ = 





(15) 

to lighten the notation. Then, given the values of /i, v, ■ • ■ at 
the vertices of a grid, their values at any point inside that 
grid can be obtained using trilinear interpolation just as we 
have done for p in Section 2.2.1. 

Now, for the configuration depicted in Fig. 1 we have, 
after some exercise of geometry, 



ViV 2 $ 



V X V X $ 



+X 2 C sm sin 2 I 



X w sin ip sin 26, 
1 



X to cos ip sin 26 



(y — /j.) sin 2ip + £ cos 2tp 



ix 2 sin 26» 

+X 2 cos 2 6 (zu sin ip — to cos ip) , 

2 2 2 2 2 

H cos ip cos 6 + v sin ip cos 6 + r\ sin 6 
+£ sin 2ip cos 2 6 + w cos ip sin 26 
+u sin ip sin 26. (16) 

Note that the above expressions are all linear in p,v,- ■ ■, 
making the situation quite simple. As an example, for 
V?$ = V^i* + V 2 V 2 ® we have 



(17) 



(18) 



2 2 / 2 2 \ 

X cos 6M^i sin ip + iycos ip — £sm2ip) 



ViVi$ = 

2 / 2 2 2 2 2 

V2V2 c 3? = x (mcos i/' sin 8 + vsin i/>sin 6 + r\ cos 



(sin 2 ip + cos 2 ip sin 2 6>) fi 

( 2 2 2 \ 2 

cos ip + sin ip sin" 6) v 4- r\ cos 6 

— £ sin 2tp cos 2 6 — uj cos tp sin 26 — uj sin tp sin 26, 

as is consistent with Castro, Heavens & Kitching (2005), 
and because ip, 6 are constants for a given ray 

Xs 

— (sin 2 ip + cos 2 ip sin 2 (fi) — (0 sin 2ip cos 2 6 
+ (cos 2 ip + sin 2 tp sin 2 (v) + (r/) cos 2 6 
— (vj) cos tp sin 26 — (oj) sin ip sin 26, 

where 

Jo Xs 

(and similarly (is), ■ ■ ■) are computed exactly as in Eq. (12). 
Note that we only need to compute {(j), ■ ■ ■ during the iV- 
body simulations and multiply appropriate coefficients as in 
Eq. (17) to obtain k finally. The components of the shear 
field (71,72) could be computed using the same formula 
as Eq. (17), but with V|$ replaced with V\<E> - V 2 2 $ 
and V 1 2 $ correspondingly using the expressions given in 
Eq. (16). 



3 iV-BODY AND RAY-TRACING 
SIMULATIONS 

To test our algorithm, we have performed a series of iV-body 
simulations for a concordance cosmology using the publicly 
available code MLAPM (Knebe, Green & Binney 2001). As it 
is not our intention to carry out very high-resolution simu- 
lations here, we only use the particle-mesh part of MLAPM so 
that our simulation grid is not self-adaptively refined. We 
have also developed a C code, RATANA (which stands for AN- 
Alytic RAy- Tracing) , to compute the convergence and shear 
fields on-the-fly as described in the above section. This sec- 
tion is devoted to a summary of our results. 

3.1 Specifications for A-body Simulations 

We consider a concordance cosmology with cosmolog- 
ical parameters Q m — 0.257, Qa = 0.743, h = 
H /(100 km/s/Mpc) = 0.719, n s = 0.963 and cr 8 = 0.769. 
The simulations start at an initial redshift z; = 49.0, and 
initial conditions (i.e., initial displacements and velocities of 
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Figure 3. Convergence and shear maps (5° X 5°) for one realisation from the tiling solutions. The convergence field re is shown as a 
colour-scale plot and the values are indicated by the colour bar below; the shear field (71,72) is shown as a flow plot and superposed on 
the k field for comparison. As expected, the shear field is tangential around high-re regions. Note that: (1) the re field in the left panel 
is computed according to Eq. (14) with the last three correction terms incorporated as described in Section 2.2.1; (2) the re field in the 
right panel is computed using Eq. (17); (3) the (71,72) field in the left panel is computed using Eq. (17) but with V^4> replaced by 
V?3? — V 2 ^ (for 71) and Vi V2$ (for 72); (4) the (71, 72) field in the right panel is computed indirectly by Fourier transforming re(6>) to 

/;2_,2 \ 

re(l), computing (71,72) = I A . ,js re, a * ^ re and finally inverse Fourier transforming (71,72) to (71,72)- 

V 1+2 ! l+ ( 2 / 



100.0 



' CO 




0.1 1.0 10.0 

k (h/Mpc) 



Figure 2. Plotted is the A^(fc) = k 3 P(k)/ (2tt 2 ), in which P(k) 
is the matter power spectrum, as a function of the wavenumber k 
in units of h Mpc -1 . The symbols with error bars represent aver- 
aged results at z = from 10 realisations for the B = 80k" 1 Mpc 
simulations. The solid curve is the corresponding result using the 
Smith et al. (2003) fit and the same set of cosmological parame- 
ters. 



particles) are generated using GRAFIC (Bertschinger 1995). 
In this work we only consider a source redshift z 3 = 1.0, 
though other values of z s or even multiple source redshifts 



can easily be implemented. The field-of-view is 5° x 5°, and 
we trace 1024 2 light rays. 

A source at redshift z a = 1.0 is about 2374/i _1 Mpc 
away from us (z = 0) in terms of comoving angular diam- 
eter distance, and it is unrealistic for us to have a simu- 
lation box which is large enough to cover the whole light- 
cone. In this work we adopt the tiling scheme introduced by 
White & Hu (2000). They use multiple simulation boxes to 
cover the light-cone between z = and z 3 , and the sizes of 
the simulation boxes are adjusted so that smaller boxes are 
used as the light rays get closer to the observer. It has been 
argued that the use of multiple tiling boxes can compensate 
the lack of statistical independence of fluctuations caused by 
using the same simulation box repeatedly. Also the variable 
box sizes mean that one can get better angular resolutions 
by using smaller boxes near the observer. 

Similar to White & Hu (2000), we choose six different 
box-sizes and 20 tiles between z — and z s , and the details 
are summarised in Table 3.1. For the A-body simulations 
(regardless of the box sizes), we use a regular mesh with 
512 x 512 x 512 cubic cells. We use the triangular-shaped 
cloud (TSC) scheme to assign the matter densities in the 
grid cell, and to interpolate the forces (Hockney & Eastwood 
1981; Knebe, Green & Binney 2001). Given the matter den- 
sities in the cells, the gravitational potential $ is computed 
using fast Fourier transform (FFT), and the gravitational 
forces (first derivatives of $) as well as the second deriva- 
tives of $ are then obtained by performing finite differences. 
These derivatives of $ are subsequently utilised by RATANA 
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Table 1. The tiling solution of our ./V-body simulations. Here 
dout is the scale factor at the time when the light rays which are 
traced leave a given tile, and B is the size of the simulation box 
in units of h^ 1 Mpc. Each simulation uses exactly 400 time steps 
from z = 49 to z = 0. AT rea j is the number of realisations for each 
value of box size. To obtain a tiling solution we randomly pick out 
two different simulation boxes with B = 240, two with B = 200, 
two with B = 160, two with B = 120, two with B = 100 and 10 
with B = 80 - a total of 20 simulation boxes of different sizes. 



(lout B(h 1 Mpc) -/V rca i ftout Bih- 1 Mpc) 7V real 
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10 
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80 
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80 
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80 
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80 
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80 




0.711 
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80 
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0.951 


80 




0.757 
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10 


0.976 


80 




0.780 
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Multipole t 

Figure 4. The convergence power spectrum A 2 . = l 2 Ci/(2n) 
as measured from our ray-tracing simulations (symbols with er- 
ror bars). The result is obtained by averaging 120 realisations 
of the tiling solution. The solid curve is again obtained using the 
Smith et al. (2003) fit of the matter power spectrum (Kaiser 1992; 
Jain & Scljak 1997), and the filled band illustrates the expected 
observational uncertainty from DES. 

to compute the convergence and shear fields as described in 
the above section. 

Note that unlike in many other works, we use the same 
grid for both iV-body and ray-tracing simulations. The TSC 
scheme we are using then results in some small-scale details 
of the matter distribution being smoothed out, as compared 
to the conventional nearest grid point (NGP) or cloud-in-cell 
(CIC) density-assignment schemes We will comment on this 
point later. 

5 In the TSC scheme, the density on a grid cell depends on the 
distribution of particles on all the 26 neighbouring grid cells; in 
the CIC (NGP) scheme, it depends on the matter distribution on the 

6 direct neighbouring grid cells (the particles in that cell only). 



3.2 Numerical Results 

In this subsection we summarise the numerical results from 
our Af-body and ray-tracing simulations. 

In Fig. 2 we compare the matter power spectra (or 
equivalently A^(fc) defined in the figure caption) computed 
from our ./V-body simulations (box size 80ft" 1 Mpc) to the 
prediction of the analytic fitting formula of Smith et al. 
(2003). We can see a good agreement, except in the range 
of 0.5h Mpc -1 < k < 2.0/i Mpc" 1 where the iV-body simu- 
lations predict a slightly higher power. However, the agree- 
ment becomes poor for k > 5.0/i Mpc -1 because the reso- 
lution of our simulations is not high enough, but this could 
be overcome in future higher-resolution simulations. 

To show that our ray-tracing simulations produce rea- 
sonable results, we first consider the convergence and shear 
maps from a chosen realisation of tiling solution, and these 
are shown in Fig. 3. We have computed the convergence 
field k(£), using the two methods outlined in Sections 2.2.1 
(Method A, Eq. (14), left panel of Fig. 3) and 2.2.2 (Method 
B, Eq. (17), right panel of Fig. 3). The two methods give 
almost identical results, and as we have checked, the differ- 
ence is in general well within the percent level. The shear 
field (71,72) is also calculated using two methods: method 
A using an equivalence of Eq. (17) as described in the fig- 
ure caption (left panel), and method B which is often used 
in the literature, namely by Fourier transforms of the con- 
vergence field which only works in the weak-lensing regime 
(right panel). The shear fields are shown in rods along with 
the convergence map shown as images. Again, the agree- 
ment for the shear field is very good, indicating that our 
ray-tracing algorithm works well. 

We also show the lensing convergence power spectrum 
measured from our ray-tracing simulations in Fig. 4. Due 
to our limit field of view of (5° x 5°), we cannot measure 
the spectrum at multiple moment £ < 100. Also, there is a 
rolloff of power at £ > 2000, which is because (1) the resolu- 
tion for our Af-body simulations is not high enough and (2) 
the TSC density- assignment scheme smooths out the small- 
scale structure more than the CIC and NGP schemes do. Both 
factors tend to suppress the convergence power spectrum at 
high t and we hope to solve this problem by using higher res- 
olution simulations and more suitable interpolation schemes, 
which is left for our future study. Otherwise, we find that 
the ray-tracing result agrees reasonably well with the ana- 
lytic prediction using the fitting formula for matter power 
spectrum by Smith et al. (2003) in some £ range, i.e., 100 < 
£ < 2000. On some scales, we see that the numerical result is 
slightly higher than the theoretical prediction. This is how- 
ever as expected because we have seen from Fig. 2 that the 
A^-body simulations give a higher matter power spectrum 
than the Smith et al. (2003) fit on some scales. The fact that 
a difference in the matter power spectra from simulations 
and analytic fitting could cause differences in the computed 
convergence power spectra has been reported and discussed 
by many authors, e.g., Vale & White (2003); Hilbert et al. 
(2009); Pielorz et al. (2010). In Fig. 4, we overplot the ex- 
pected observational uncertainty from DES using the survey 
parameters / s k y ~ 0.12, n g = 10/arcmin 2 , 7 ; 2 nt ~ 0.16 where 
/sk y ,n g and 7 ; 2 nt denote the sky coverage, number of galax- 
ies per arc-minute squared and the mean-square intrinsic 
ellipticity, respectively. 
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Note that in our numerical simulations we have not in- 
cluded the lens-lens coupling and second-order corrections 
to the Born approximation. In Appendix A we will outline 
how these can be incorporated in future higher-resolution 
simulations. 



4 DISCUSSION AND CONCLUSION 

The correlations in the distorted images of distant galaxies, 
induced by cosmic shear, hold information about the distri- 
bution of matter on a wide range of scales in the universe. 
In order to take full advantage of current and future weak 
lensing data sets to constrain cosmology, using information 
from both the linear and non-linear regimes, one needs a so- 
phisticated algorithm to measure the shear and convergence 
fields from A-body simulations, and to construct statistical 
quantities. This is traditionally done using the 'plane-by- 
plane' discrete lens-plane algorithm - trace the virtual light 
rays and calculate the deflection caused by the density field 
projected onto a number of 2-D lensing planes. 

In this work, we propose an improved ray-tracing algo- 
rithm. We calculate the deflection of the light rays caused 
by the detailed 3-D density fields living on the natural sim- 
ulation mesh, rather than the simplified density distribution 
projected onto some 2-D planes. We evaluate the shear and 
convergence fields by analytically integrating the deflection 
as the light rays go through the individual simulation grid 
cells. This approach is easy to implement and computation- 
ally inexpensive. It avoids numerical integration, and expen- 
sive data storage since it is performed on the fly. We apply 
the algorithm to our simulations, and find good agreement 
with the Smith et al. (2003) fit, and consistency with the 
published results in Sato et al. (2009). 

The on-the-fly l.o.s. integration is computationally eco- 
nomic. In the RATANA code, most computation time is spent 
on the A-body part. Suppose Aj is the number of grid cells 
in our mesh, then the FFT requires 3A| log 2 Nd operations 
each time step, not including other operations such as differ- 
encing the potential to obtain the force on the mesh, assign- 
ing particles and computing densities on all the grid cells 
and particle movements. In contrast, if we let Ai OB = Nd 
(which is enough for accuracy), then there are only Aj rays 
to trace, and for each ray we have ^ 10 2 operations. We 
have checked the simulation log file and found that there is 
no significant difference in the times used by each step before 
and after the ray-tracing part of RATANA has been triggered. 

Analytic formulae are often more useful than purely nu- 
merical results in tracing the physical contents of a theory. 
For example, in Eqs. (12, 13), it is easy to check which 
terms contribute the most to the final result: obviously, 
in the small-angle limit, i.e., 0,tp <C 1, terms involving 
dijdi, and a large part of d^ could be neglected because 
sin 6, sin ip <C 1; also at least one of a,b,c vanishes and 
abc = for all grid cells, further simplifying di,d2\ further- 
more, terms in Eq. (12) with coefficient 1/xs contribute lit- 
tle because Xs Xu,i- Such observations can be helpful in 
determining which terms have important effects in certain 
regimes. 

Note that the dependence on Xs [cf. Eq. (12)] could be 
taken out of the analytical integration, meaning that the al- 
gorithm can be straightforwardly generalised to include mul- 



tiple source redshifts with very little extra computational 
effort (mainly in determining where to start the integration 
for a given source redshift). The algorithm can also be easily 
generalised to compute the flexion, which depends on higher- 
order derivatives of the lensing potential, and is expected to 
give more accurate results than the multiple-lens-plane ap- 
proximation. 

The algorithm has many other flexibilities too. As an 
example, the analytic integration of the projected den- 
sity and potential fields along the l.o.s. can be performed 
on an adaptive rather than a regular grid with careful 
programming, which means that higher resolution can be 
achieved in high density regions, as in the adaptive PM sim- 
ulations. Also, the analytic integration can be easily gener- 
alised to other algorithms to compute the 3-D shear field 
(Couchman, Barber & Thomas 1999). 

We also give prescriptions to include second-order cor- 
rections to the results, such as the lens-lens coupling and 
corrections to the Born approximation, in Appendix A. It is 
interesting to note that, by running the A-body simulations 
backwards in time, we can still compute the convergence and 
shear fields on-the-fly even if the light rays are not straight. 

To conclude, the algorithm described here is efficient 
and accurate, and is suitable for the future ray-tracing sim- 
ulations using very large A-body simulations. It will be in- 
teresting to apply it to study the higher-order statistics of 
the shear field and the lensing excursion angles, and these 
will be left for future work. 
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APPENDIX A: BEYOND THE FIRST-ORDER APPROXIMATIONS 

In the attempt to trace light rays on the fly, we set up a bundle of l.o.s. before the JV-body simulation starts. But because we 
do not know the exact paths of those light rays which finally end up at the observer, we have to assume that they are straight 
lines even though they are not in reality. This so-called Born approximation is generally quite good in the weak lensing regime, 
but can lead to non- negligible errors on small scales (Hilbert et al. 2009). Furthermore, in the above treatment we have also 
neglected the lens-lens coupling, which accounts for the fact that the lenses themselves (the large-scale structure) are distorted 
by the lower-redshift matter distribution. 

Hilbert et al. (2009) take account of the lens-lens coupling and corrections to the Born approximation using the multiple- 
lens-plane approximation. In such an approach, the light rays get deflected and their paths are recomputed when and only 
when they pass by a discrete lens plane. 

Since our algorithm goes beyond the discrete lens-plane approximation and is able to trace the detailed matter distribution, 
we want to generalise it to include those corrections as well. In this Appendix we shall derive an analytical formula for the 
distortion matrix with the lens-lens coupling taken into account, and describe how the corrections to the Born approximation 
can be incorporated as well. 

Obviously, to go beyond the Born approximation, the light rays are no longer straight and thus the l.o.s cannot be set 
up before the iV-body simulation has finished. Instead, we have to start from the observer today and go backwards in time to 
compute the distortion matrix Eq. (3). We shall discuss below how this could be realised in practice, but at this moment let 
us simply assume that we can go backwards in time, and know the value of the lensing potential "I> and its derivatives along 
the l.o.s.. 



Al Corrections to the Born Approximation 

The corrections to the Born approximation are easy to implement. According to Eq. (2), the total deflection of a light ray 
is the sum of the deflections by the matter in each grid that ray passes on its way towards the lensing source. Suppose |^™^ 
denotes the value of £ after the light ray crosses the n-th grid on its way (n increases with the distance from the observer, 
n — 1 corresponds to the grid which the observer is in, and £} 0) = Co), then 

rxi n) 

& } = ^ n " 13 - 2 / — — -V,-$dx, (Al) 
J X M XXs ? 

where xi™' = rnin |xu*> Xs^ | an d Xi™' = max |x; s , Xa*' m which Xi? > xT are respectively the \- values at the two ends of 

the current time step, and Xb > Xa the x-values of the two intersections between the light ray and the n-th grid. Using the 
expressions given in Sect. 2.2.2, it is easy to write V^"!? and V^ 2 $ in terms of polynomials of \. Then the above integral can 
be performed analytically as before. In this way, each time the light ray crosses a grid, we update its orientation according to 
the above equation, and thus the corrections to the Born approximation can be incorporated. 

Note that in this approach the light rays are deflected many more times than in the multiple-lens-plane approximation 
and the detailed matter distribution has been fully taken account of. 



A2 Lens-lens Coupling 

As mentioned earlier, the lens-lens coupling has been neglected in the above treatment because in Eq. (3) we have used the 
approximation V xoi ~ V Xi . Let us now have a look at what happens when this approximation is dropped. 
Note that in the expression 

<S ^%r^V foi V<y<i> (x,£) d x + 6a, (A2) 
x 

the argument of $ is £ while one of the derivatives is with respect to £o- We can utilise the chain rule to write V^ oi = 
(Vf oi £j) V$. = AijV '| . where we have used the definition of Ay given in Sect. 2.1. Then the above equation becomes 

A\ (x„|) = S'j - 2 / 5 (x,Xs)V l V fc $ ( X ,|) A k j ( X ,£) d X , (A3) 
Jo 

where for simplicity we have used V; = Vj 4 . With the A k j term in the integrand, Eq. (A3) now includes the lens-lens coupling, 
and will be our starting point here. 

Again, let us consider the integral in Eq. (A3) after the light ray crosses the n-th grid on its way towards the lensing 
source. The discrete version of Eq. (A3) is 

, x <,"> , . 
("O^V _ (™- 1 ) J 4 4 ._2 / X (Xs — X) (n)^fe 



j j 



M fc ,V'V fc *dx (A4) 

Xs 
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where ^ A % j is the value of A 1 j after the light ray has crossed the n-th grid, and (°'A % j = S l j as is easy to see. This formula 
has three advantages as compared to the multiple-lens-plane approximation: 

(i) As before, the light rays between z = and z E are divided into many more segments, and the fine structure of the 
matter distribution is included naturally, without squeezing the matter and using impulse approximations. 

(ii) As will be shown below, the integration can be evaluated analytically rather than numerically. 

(iii) Note that we can use ( n 'A k j rather than ( n ~ 1 > A k j in the integrand, which will give more accurate results, because 
using ( n ~ 1 ^A k j would mean that the contribution to the lens-lens coupling from the matter in the n-th grid is ignored. In the 
multiple-lens-plane approximation which typically uses 20 ~ 30 lens planes, the n-th plane could contain a significant amount 
of matter, and neglecting its contribution could make the results less accurate. 

Eq. (A4) is exact, but we only want the result to second order in VV<3>. Therefore we can iterate once and write an 
approximate solution as 
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Following the approach taken in Sect. 2.2.1 we can write 
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(A6) 



where Xa ls defined in Eq. (Al), and Kjv (N £ {1, 2, 3, 4}) is a 2 x 2 matrix whose ij-component depends on the orientation 
of the l.o.s. segment inside the n-th grid (where it is taken to be straight) and the values of V l Vj<3? at the vertices of the n-th 
grid. Note however that Kjv is independent of x- The expressions are similar to the djvs defined in Sect. 2.2.1 and we shall 
not write them explicitly here. 

Substituting Eq. (A6) into Eq. (A5), we find 
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in which we have written (again, by defining x = X~ Xa , X = X ~ X ( a and Xu} = xl™ - Xa) 
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The above expressions look rather heavy, however, they are analytic and as a result are very easy to implement in the 
ray-tracing simulation codes, by writing functions that take M, N, Xa , Xu > X; as parameters and return I\, I2 as outputs. 
Furthermore, since the grid size (< 0.2/i _1 Mpc) in the iV-body simulations is small enough compared with the typical inter- 
plane distances in the multiple-lens-plane approximations (10 ~ lOO/i -1 Mpc), we can drop the I2(N,M) terms to a very 
good approximation, which will greatly simplify the results. 

Note that the distortion matrix A Z j computed in this way is not symmetric because of the matrix multiplications. However, 



using Eq. (5), it is straightforward to compute 72 = 
u = (A\ - A\) /2 



[A 1 2 + A 2 ij jl. In addition, we could also calculate the rotation u as 



A3 Going Back In Time 

As mentioned above, to include the actual deflections of the light rays which end up at the observer, we have to start from 
the observer and go backwards in time until encountering the source. This obviously can only be done after the iV-body 
simulation has finished. 

One way to go backwards is to record the information about the gravitational potential $ and its derivatives in a light 
cone during the simulation, and then post-process the light-cone data. This means that a large amount of dump data has to 
be stored. 

Alternatively, one can think of running the iV-body simulation "backwards". To be more explicit, the simulation is first 
run in the forward direction from a high redshift until today, and we obtain the particle positions and velocities at present; 
then we reverse the directions of the gravitational force and the particle velocities, and evolve the system back until z s using 
the same time-stepping scheme as in the forward simulation. In this way, the actual light rays and distortion matrix could be 
built up on the fly, and there is no need to store a lot of dump data. 



